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SECTION  1 


INTRODUCTION 

This  report  describes  the  application  of  a  large-scale  numerical  wave 
solver  to  the  time  domain  study  of  seismic  wave  phenomena  in  filled  basins. 

The  solver  uses  an  explicit  finite  element  algorithm  designed  for  fully  vec¬ 
torized  execution  on  the  CRAY-1,  permitting  the  solution  of  discrete  hyperbolic 
problems  on  a  scale  at  least  one  order  of  magnitude  faster  than  conventional 
scalar  machines  allow  (e.g.,  CDC  7600).  This  computing  power  is  sufficient 
to  analyze  large  2-D  inhomogeneous  sections  of  the  Earth's  upper  crust  with, 
for  example,  a  minimum  of  2  Hz  frequency  resolution  of  Rayleigh  and  shear 
waves  in  a  50  km  x  7  km  grid  (120  x  1000  elements).  Such  geologic  models  are 
essential  to  the  present  study  of  body  and  surface  wave  propagation  in  realis¬ 
tic  basin  structure,  particularly  the  questions  of  edge  effects,  Rayleigh  wave 
interactions,  and  wave  coupling  between  basins  on  either  side  of  a  mountain. 
Because  large  scale  computing  is  fairly  new,  the  numerical  modeling  aspects 
of  the  problem  are  described  in  detail  with  emphasis  on  wave  solver  mechanics. 

In  Section  2,  a  set  of  four  basin  configurations  are  analyzed  for  two  load¬ 
ing  conditions  and  two  velocity  functions.  The  three  single  basin  models 
include  a  steep  faulted  flank,  an  echelon  faulted  flank  and  a  dipping  flank.  A 
larger  scale  basin-range-basin  model  completes  the  set.  Full-field  synthetic 
seismograms  are  displayed  and  interpreted  for  each  model.  In  Section  3,  the 
finite  element  theory  is  described  for  a  cartesian  computational  mesh.  The 
emphasis  is  on  methods  of  reducing  the  operations  count  to  minimize  computer 
processing  time.  In  Section  4,  the  calculations  and  theory  are  summarized  and 
conclusions  drawn. 
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SECTIOM  2 


FULL-FIELD  SEISMIC  RESPONSE  IN  BASINS 

The  two-dimensional  numerical  calculations  described  in  this  section 
include  a  layered  halfspace  model  and  three  symmetric  basin  models  indicative 
of  local  geologies  in  the  Basin  and  Range  province  of  the  western  United 
States.  In  addition  a  larger  model  of  two  basins  separated  by  a  mountain 
range  is  used  to  investigate  the  degree  of  seismic  coupling  between  basins. 
Near-surface  seismic  velocity  functions  are  constructed  using  typical  data 
compiled  by  Battis  (1981)  from  a  variety  of  Basin  and  Range  studies.  Velocity 
gradients  in  bedrock  are  taken  from  Prodehl  (1979),  based  on  USGS  crustal 
refraction  data.  Vertical  and  horizontal  surface  loadings  on  the  model 
centerline  provide  the  seismic  input,  while  vertical  velocity  seismograms 
and  particle  motion  plots  are  used  to  interpret  response. 

2. 1  Geologic  Models 

The  Basin  and  Range  province  is  characterized  by  east-west  lateral 
extension,  crustal  thinning  and  uplift  forming  many  parallel  faults 
trending  north-south.  The  resulting  topography  consists  of  parallel,  elon¬ 
gated  basins  separated  by  narrow  mountain  ranges.  The  formative  mechanism 
is  thought  to  be  block  subsidence  due  to  extension,  Stewart  (1971),  although 
secondary  block  tilting  may  also  be  significant.  A  schematic  illustration  is 
given  in  Fig.  2-1.  The  resulting  horst  and  graben  model  involves  high  angle 
faulting  with  typical  dips  between  45  and  70  degrees  and  fault  throws  on  the 
order  of  a  kilometer  or  more.  Average  basin  depth  is  2  km,  width  is  15  to 
20  km,  and  length  is  approximately  70  km.  The  intervening  mountain  ranges 
are  usually  about  15  km  wide  and  maximum  elevation  above  basin  fill  may  range 
from  0.3  to  1.2  km.  These  typical  data  were  compiled  by  Battis  (1981)  in  a 
study  of  the  open  literature  on  Basin  and  Range  geology  and  tectonics. 
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Figure  2-1. 


A  schematic  of  the  graben  structure 
caused  by  deep  seated  extension  and  block 
subsidence  in  the  Basin  and  Range  province. 
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Figure  2-2.  Nominal  Basin  and  Range  graben  basin,  16  km  wide 
and  2  km  deep  with  faulted  flanks  dipping  at 
approximately  60  degrees. 
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Based  on  the  above  data,  a  nominal  basin  is  assumed  to  be  2  km  deep  and 


16  km  wide,  with  normal  faulted  flanks  dipping  at  approximately  60  degrees. 

The  length  is  assumed  to  be  four  times  the  width — sufficient  to  justify  a  two- 
dimensional  model  of  the  structure  for  numerical  analysis.  Such  a  basin  is 
illustrated  in  Fig.  2-2  for  the  asymmetric  faulting  case.  This  cross  section 
suggests  the  models  for  numerical  study  shown  in  Fig.  2-3.  The  simple  layer 
model.  Fig.  2- 3a,  establishes  a  baseline  for  response  in  the  basin  fill,  with¬ 
out  edge  effects.  The  basin  models.  Fig.  2-3b,c,d,  provide  three  types  of 
edges — a  single  faulted  flank  (b),  an  echelon  faulted  flank  (c),  and  a 
shallow  dipping  flank  (d)  which  is  the  limiting  case  of  echelon  faults.  The  0 
models  assume  symmetry  across  the  basin  axis  in  order  to  minimize  problem  £  i 
Some  effects  of  asymmetry  can  be  estimated  from  them.  On  a  larger  scale, 

Fig.  2-4  shows  the  two-basin  model  with  an  intervening  mountain  range.  Her 
and  width  of  the  range  are  taken  as  1  km  and  16  km,  respectively. 

2. 2  Velocity  Models 

The  graben  basins  have  been  continually  filled  during  their  development — 
predominantly  with  alluvium  from  stream  erosion  of  neighboring  mountains,  but 
with  significant  ash  and  lava  deposits  during  periods  of  volcanism.  The  sur¬ 
rounding  crustal  blocks  include  a  variety  of  sedimentary,  volcanic,  and  plu- 
tonic  rock  types.  This  situation  indicates  a  highly  inhomogeneous  velocity 
distribution  for  the  basin  fill  and  bedrock. 

Basin  seismic  velocities  are  tabulated  by  Battis  (1981)  from  studies  by 
Thompson,  et  al.  (1967)  and  Healv  and  Press  (1964).  The  average  P-wave  veloc¬ 
ities  and  bounds  are  shown  in  Fig.  2- 5a.  Superposed  on  these  is  a  linear 
approximation  to  the  average,  starting  at  2  km/sec  on  the  free  surface  with 
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a.  LAYER  b.  FAULTED  FLANK 


c.  ECHELON  FAULTED  FLANK 


d.  DIPPING  FLANK 


Figure  2-3.  Layer  and  basin  cross  sections  for  numerical  analy¬ 
sis  of  edge  effects.  Basins  are  2  km  deep  and  8  km 
from  centerline  to  edge. 


BASIN  -RANGE-BASIN 


Figure  2- A. 


Two-basin  model  with  an  intervening  mountain  range, 
16  km  wide  and  1  km  high. 
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BASIN  VELOCITY  MODEL  b.  UPPER  CRUST  VELOCITY  MODEL 


Figure  2-5.  Velocity  functions  assumed  for  basin  fills  and 
surrounding  bedrock.  Basin  data  compiled  by 
Battis  (1981)  and  crustal  gradients  derived 
from  Prodehl  (1979). 


a  1  km/sec/km  gradient  down  to  1  km  depth  and  a  0.5  km/sec/km  gradient  to 
2  km  where  the  velocity  is  3.5  km/sec.  This  approximation  ignores  the  low 
speed  surface  layer  with  velocity  between  1  and  2  km/sec  in  the  upper  0.2  km. 
Of  course, the  average  data  does  not  indicate  lateral  variations  in  velocity 
across  the  basin  due  to  differences  between  interior  (flats  or  playa)  and 
perimeter  (alluvial  fan)  deposits;  nor  does  it  include  volcanic  anomalies 
like  intrusions,  old  surface  flows,  etc.  These  will  be  neglected  here  in 
comparison  to  the  lateral  change  between  alluvium  and  bedrock  properties  at 
the  basin  edge. 

Bedrock  seismic  velocities  have  not  been  studied  to  the  extent  that 
basin  velocities  have;  however,  estimates  from  in-situ  measurements  of 
exposed  bedrock  indicate  a  range  from  4.7  to  5.1  km/sec,  Thompson  et  al. 
(1967).  The  variation  of  bedrock  velocity  with  depth  is  not  known.  USGS 
crustal  refraction  data  in  the  Basin  and  Range,  reduced  by  Prodehl  (1979) 
yield  a  crude  estimate  of  the  upper  crust  gradient  at  0.4  km/sec. km.  This 
is  based  on  a  surface  velocity  of  4  km/sec  varying  smoothly  to  6  km/sec  at 
5  km  depth,  from  Prodehl's  velocity  function  drawings.  The  actual  gradient 
is  probably  not  constant  and  could  easily  be  half  the  above  value,  assuming 
a  surface  velocity  of  5  km/sec  for  example.  The  resulting  upper  crust  veloc¬ 
ity  model  is  shown  in  Fig.  2-5b.  At  2  km  depth  the  P-wave  velocity  is 
4.8  km/sec  rising  linearly  to  6  km/sec  at  5  km  depth.  Below  5  km  the  data 
show  a  very  low  gradient,  hence  a  constant  velocity  is  assumed.  A  piecewise 
constant  model  is  also  drawn  in  the  figure,  with  2.8  km/sec  in  the  basin  fill 
and  5  km/sec  in  the  rock.  This  choice  of  basin  velocity  gives  the  same 
vertical  transit  time  as  the  piecewise  linear  model. 

The  remaining  seismic  data  needed  for  analysis  are  S-wave  velocities 
and  material  densities.  A  good  S-wave  approximation  for  the  deeper  basin 
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fill  and  bedrock  is  the  Poisson  assumption,  v  =  0.25  giving  a  ratio  of  P  to 


S  velocity  of  Vp/Vg  =  *3.  For  simplicity,  this  will  also  be  used  for  the 

shallow  basin  fill,  although  data  show  that  the  ratio  is  generally  higher 

3 

for  the  weathered  layer.  Densities  are  assumed  to  be  2.2  gm/ cm  for  alluvium 

3 

and  2.5  gm/cm'  for  bedrock.  These  values  are  probably  within  10  to  15  percent 
of  the  actual  depth  dependent  densities. 


2 . 3  Finite  Element  Models 

The  geologic  cross  sections  in  Figs.  2-3,4  are  modeled  by  cartesian, 
plane  strain  finite  element  grids.  A  7  x  20  km  model  for  the  echelon  faulted 
flank  case,  Fig.  2-3c,  is  illustrated  in  Fig.  2-6,  for  example.  Boundary  con¬ 
ditions  are  such  that  the  surface  is  traction  free,  the  left  side  is  a  line 
of  symmetry,  while  the  right  and  bottom  sides  are  transmitting  boundaries 
(normal  impedance  type).  Only  seismic  line  sources  on  the  surface  are  consid¬ 
ered,  consisting  of  normal  or  tangential  tractions  in  a  Gaussian  distribution 
over  a  few  nodes  about  the  model  centerline.  This  spread  of  surface  loading 
is  necessary  to  minimize  nonphysical  element  deformational  modes  associated 
with  single  node  forcing  functions. 

Element  size  depends  on  the  miminum  wavelength  to  be  resolved.  From  the 

velocity  functions  described  above,  clearly  the  slowest  phase  in  these  models 

will  be  the  Rayleigh  surface  wave,  for  which  V  -  0.919  V  >  1.06  km/sec 

K  b  — 

(based  on  V  =  V  //3  -  1616  km/sec  at  the  free  surface).  Assuming  the  maximum 

o  p 

frequency  of  interest  is  2  Hz  (1/2  second  period)  then  the  minimum  wavelength 
is  503  meters.  The  minimum  number  of  elements  per  wavelength  is  no  less 
than  10,  whence  the  maximum  element  size  is  50  meters.  To  gage  whether  10 
elements  per  wavelength  is  adequate,  it  is  sufficient  to  compare  maximum  input 
frequency  to  the  minimum  natural  frequency  of  50  meter  square  elements.  From 
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Section  3.5,  f  =  1/T  -  V  /(8.89£x)  ~  4.45  Hz,  hence  the  input  of  2  Hz  is 

less  than  half  the  element  ringing  frequency  and  adequate  for  modeling  purposes. 

Because  the  normal  impedance  type  transmitting  boundary  conditions  used 
here  are  not  perfect,  some  energy  will  be  reflected  back  into  the  grid.  To 
minimize  the  possibility  of  confusing  reflected  phases  from  the  boundary  with 
refracted  phases,  elements  are  graded  in  size  near  the  boundary  to  introduce 
dispersion  and  scattering.  Therefore,  elements  are  made  square  in  a  3  x  16  km 
strip  including  the  basin  in  Fig.  2-6,  but  are  graded  by  3  percent  per  element 
to  the  bottom  and  right  boundaries  yielding  a  maximum  element  aspect  ratio 
of  3  there.  This  is  a  crude  solution  at  best  to  the  transmitting  boundary 
problem  and  will  be  evaluated  later. 

The  resulting  finite  element  models  to  be  solved  on  the  CRAY-1  contain 
36,000  elements  for  each  of  the  basin  flank  cases  shown  in  Fig.  2-3,  and 
120,000  elements  for  the  basin-range-basin  model  in  Fig.  2-4.  The  explicit 
finite  element  code  used  (for  which  the  theory  is  described  in  Section  3)  re¬ 
quires  two  microseconds  (2  x  10  ^  sec)  to  evolve  each  element  one  timestep. 

The  timestep  is  taken  as  80  percent  of  the  minimum  element  transit  time  which, 

-3 

from  Fig.  2-5,  gives  6.7  milliseconds  (6.7  x  10  sec)  for  50  m  elements.  Finite 

element  run  times  are  chosen  to  propagate  the  free  surface  Rayleigh  wave  across 

the  model  (approximately)  requiring  1500  timesteps  for  the  constant  velocity 

flanks,  2000  steps  for  the  variable  velocity  flanks  and  5000  steps  for  the 

basin- range-basin  case.  Multiplying  number  of  model  elements  by  number  of 

6  6  8 

timesteps  gives  54  x  10  ,  72  x  10  and  6  x  10  element-timesteps,  respectively, 
for  these  cases,  requiring  108,  145  and  1200  seconds  for  their  execution. 

2. 4  Synthetic  Seismograms 

Model  response  to  the  surface  traction  inputs  is  recorded  at  many  output 
points  over  the  free  surface,  and  also  through  depth  at  selected  ranges. 
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Only  surface  response  is  examined  here,  however.  Because  the  model  is  linear 
(small  displacements),  it  suffices  to  use  a  Green's  function  formulation,  i.e., 
the  original  source  time  function  is  a  generalized  function  (delta  function, 
step,  etc.).  Any  desired  time  history  may  then  be  obtained  by  convolution 
(Duhamel  superposition).  Of  course,  frequency  resolution  is  no  better  than 
the  original  finite  element  grid  will  allow.  Thus,  from  a  single  finite  ele¬ 
ment  run,  synthetics  are  constructed  for  any  source  time  function.  The  effect 
of  instrument  response,  if  necessary,  is  included  by  a  second  convolution.  A 
change  in  source  type  or  addition  of  new  structure  to  the  model  requires  a  new 
finite  element  run. 

Before  examining  model  response,  some  comments  on  interpretation  are  in 
order.  The  suite  of  synthetics  typically  shows  a  number  of  phases  identified 
by  their  correlation  in  space  and  time  from  one  trace  to  the  next.  There 
are  a  variety  of  tools  available  to  aid  in  extracting  quantitative  information 
from  the  synthetics.  One  is  a  correlation  function  to  track  phases  (station¬ 
ary  phase  points)  automatically  yielding  travel  time  curves,  phase  velocities, 
etc;  and  another  is  beam-forming  to  determine  the  direction  of  incident  waves. 
In  the  following,  only  qualitative  results  (visual  correlations)  are  discussed 
but  it  will  be  clear  that  quantitative  tools  in  conjunction  with  finite  ele¬ 
ment  solutions  offer  major  advantages.  Another  comment  concerns  the  effect 
of  structural  modifications  in  the  model.  These  are  easily  found  by  subtract¬ 
ing  the  original  and  modified  suite  of  synthetics,  showing  seismic  influence 
independently  of  the  complete  seismogram  (which  contains  a  variety  of  strong 
phases  masking  subtle  structural  effects).  Applications  include  effects  of 
local  changes  in  wavespeed  due  to  inclusions  and  cavities,  corner  diffractors, 
etc.  Finally,  it  should  be  recognized  that  detailed  interpretation  of  seismo¬ 
grams  requires  a  ray  tracing  capability  in  inhomogeneous  media.  Therefore,  to 
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quantify  aspects  of  the  following  visual  interpretations,  a  two-point  ray  trac¬ 
ing  analysis  would  be  useful  and  some  will  be  described  in  the  final  report. 

2.5  The  Layer/Halfspace  With  Constant  Velocity 

To  establish  a  baseline  for  subsequent  basin  analyses,  it  is  useful  to 
start  with  the  simple  layer  model  loaded  by  normal  surface  traction  at  the 
origin.  This  first  calculation  assumes  the  piecewise  constant  velocity  func¬ 
tion  in  Fig.  2- 5b.  The  model  and  vertical  velocity  seismograms  are  shown  in 
Fig.  2-7a.  The  Green's  function  is  convolved  with  a  2  Hz  wavelet  for  these 
synthetics  and  amplitude  is  normalized  by  the  Rayleigh  wave  maximum  on  the 
second  trace,  300  meters  from  the  origin.  The  dominant  phase  is  clearly  the 
Rayleigh  wave  and  to  see  the  body  waves  each  trace  is  multiplied  by  10  and 
truncated  as  shown  in  the  upper  suite.  Data  was  collected  at  a  finer  spatial 
increment  but  the  300  meter  spacing  in  the  figure  is  convenient  for  visual 
interpretation. 

Slope  of  the  Rayleigh  wave  arrivals  (time-distance  curve)  in  Fig.  2- 7a 
yields  a  phase  velocity  of  approximately  1.45  km/sec,  which  is  2.4  percent 
lower  than  the  theoretical  value  for  a  halfspace  of  the  layer  medium  (basin 
fill)  only.  Because  the  wavelength  is  0.743  km  (providing  15  elements  per 
wavelength)  the  2  km  interface  is  2.69  wavelengths  below  the  surface,  hence 
the  surface  wave  is  barely  affected  by  the  faster  halfspace.  The  2.4  percent 
discrepancy  between  exact  and  measured  phase  velocity  is  caused  by  the  dis¬ 
crete  finite  element  model — a  combination  of  discrete  data  (15  elements  per 
wavelength),  the  lumped  mass  approximation,  and  reduced  stiffness  of  element 
bending  modes. 

Body  waves  are  much  lower  in  amplitude  and  more  complicated  to  interpret 
than  surface  waves.  The  body  wave  arrivals  in  Fig.  2- 7a  are  a  summation  of 
direct  waves  in  the  layer,  reflected  waves  from  the  interface  and  head  waves. 
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Spacing  between  output  points  on  the  surface  is  300  meters. 


The  direct  ?-wave  gracing  the  surface  is  the  first  arrival  out  to  a  range  of 

3.35  km,  after  which  head  waves  overtake  it.  The  head  wave  begins  at  the 

critical  range  of  2.71  km  (travel  time  =  1.67  sec)  and  is  seen  most  distinctly 

as  the  very  low  amplitude  first  arrival  beyond  9  km.  Reflected  waves  from 
the  interface  follow  the  direct  P-wave  and  are  seen  in  the  range  from  0  to 

1.5  km  with  peak  amplitude  arriving  at  about  2  sec  (first  arrival  at 

1.43  sec.).  This  arrival  is  obscured  by  the  Rayleigh  wave  from  1.5  to  4  km 

and  reappears  at  4  km  out  to  8  km  or  so,  but  beyond  amplitude  drops  consider¬ 

ably.  The  change  in  amplitude  is  caused  by  constructive  interference  of 
direct,  reflected  and  head  waves  in  the  range  from  4  to  8  km  due  to  closely 
spaced  arrival  times.  3eycnd  8  km  the  arrivals  are  separated  in  time  and 
do  not  interfere  significantly. 

The  S-wave  phases  at  any  range  follow  the  same  travel  path  as  correspond¬ 
ing  P-waves  (because  Vp/Vs  is  constant  over  the  model),  but  have  travel  time 
,7  longer.  They  are  obscured  by  the  Rayleigh  wave  out  to  8  km  and  are  diffi¬ 
cult  to  pick  thereafter,  probably  because  of  constructive  interference  by 
?-wave  multiples.  This  interference  begins  to  form  a  ducted  Rayleigh  wave 
beyond  9  km  (starting  at  6  sec.). 

Reflections  from  the  model  boundary,  although  minimized,  are  still  in 
evidence,  particularly  near  the  origin  at  4  and  8  seconds.  These  dispersed 
wavetrains  are  principally  reflected  from  the  graded  elements  rather  than 
the  boundary  itself,  and  their  influence  downrange  is  difficult  to  assess. 
However,  based  on  the  absence  of  clear  reflected  phases  prior  tc  the  Rayleigh 
wave  arrival, the  boundary  reflections  downrange  appear  to  be  minimal.  Behind 
the  Rayleigh  wave  this  may  not  be  the  case. 

2 . 6  The  Layer/Half space  With  Linear  Velocity 

The  second  layer  calculation  uses  the  same  finite  element  model  as  above 
but  with  the  piecewise  linear  velocity  function  in  Fig.  2-5a;  again  loaded  by 
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normal  surface  traction  at  the  origin.  Normalized  seismograms  for  the  2  Hz 
wavelet  are  shown  in  Fig.  2- 7b.  The  Rayleigh  wave  is  clearly  more  dispersed 
in  this  case,  caused  in  large  part  by  the  velocity  gradient — but  also  increased 
by  coarser  frequency  resolution  near  the  free  surface  (10  elements/wavelength 
rather  than  16  as  in  Fig.  2- 7a).  The  dispersion  is  readily  seen  by  viewing 
obliquely  along  the  crest  of  the  phase  in  Fig. 2- 7b.  At  any  range,  phase 
velocity,  measured  by  the  local  slope  of  a  line  connecting  stationary  phase 
points,  is  approximately  constant;  but  time  of  peak  amplitude  arrival  is 
continuously  retarded. 

Body  waves  are  seen  to  be  stronger  phases  in  this  case  due  to  the  linear 
gradient.  The  arrivals  now  consist  of  refracted  waves  in  the  layer,  reflected 
waves  from  the  interface  discontinuity  and  refracted  waves  in  the  underlying 
bedrock.  The  piecewise  linear  velocity  gradient  curves  the  wavefronts  and 
rays  and  reduces  the  head  wave  to  insignificance  in  comparison  to  the  analo¬ 
gous  refracted  wave.  The  predominant  first  arrivals  in  the  range  out  to 
4  km  are  refracted  P-waves  in  the  layer.  Between  4  and  7.5  km  the  arrivals 
are  strengthened  by  constructive  interference  of  reflected  and  refracted 
bedrock  waves  in  addition  to  refracted  basin  waves;  and  beyond  7.5  km  the 
arrival  is  predominantly  refracted  through  bedrock  (corresponding  to  head 
waves  in  the  constant  velocity  model). 

Distinctive  S-wave  phases  are  more  difficult  to  identify  because  of 
interference  by  other  arrivals — Rayleigh  waves  close-in  and  constructively 
interfering  P-wave  multiples  further  out.  Refracted  basin  S-waves  are  clear¬ 
est  at  intermediate  range,  between  6  km  and  7,5  km  for  example,  with  4  to  5 
second  arrival  time.  Refracted  bedrock  S-wave  arrivals  continue  smoothly 
from  the  basin  phase  beyond  7.5  km  but  are  preceeded  by  a  P-wave  phase 
(between  10.5  and  13.5  km)  and  followed  by  an  S-wave  phase  (both  identified 
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on  che  basis  of  travel-time  slope,  i.e.,  phase  velocity).  These  are  probably 
the  result  of  multiple  ?-  and  S-wave  refractions  and  mode  conversion, 
although  a  detailed  ray  analysis  would  be  necessary  to  confirm  this. 

Comparing  seismograms  in  Fig.  2- 7a  and  2- 7b  close  to  the  origin  shows 
that  reflections  off  the  bottom  model  boundary  are  lower  amplitude  and  less 
coherent  for  the  linear  velocity  structure.  This  follows  because  the  element 
grading  effect  is  similar  to  that  produced  by  the  velocity  gradient.  In 
addition,  the  curving  of  rays  (energy  travel  paths)  in  the  model  turn  much 
of  the  downward  propagating  energy  back,  towards  the  surface.  Based  on  these 
observations  the  grading  and  transmitting  boundary  conditions  appear  to 
perform  adequately  for  the  linear  velocity  models. 

2. 7  Basins  With  Constant  Velocity 

Analysis  of  basin  response  will  begin  by  considering  the  piecewise  con¬ 
stant  velocity  function  in  the  three  basin  models.  The  models  and  2  Hz 
synthetics  are  displayed  in  Fig.  2~8a,b,c.  Basin  flanks  are  resolved  in 
the  stepwise  fashion  indicated  in  Fig.  2-6.  In  comparison  to  the  layer 
results,  the  distinctive  feature  of  these  seismograms  is  transmission  and 
reflection  of  the  Rayleigh  wave  at  the  basin  boundary  8  km  from  the  center- 
line  (line  of  symmetry).  Comparing  the  magnified  seismograms  shows  that  the 
steep  faulted  flank  (Fig.  2-8a)  is  the  most  effective  reflector,  the  echelon 
faulted  flank  (Fig.  2-8b)  is  somewhat  less  effective,  and  the  dipping  flank 
(Fig.  2-8c)  is  a  relatively  poor  reflector.  However,  in  terms  of  transmission, 
the  steep  and  dipping  flanks  are  comparable,  and  both  are  better  transmitters 
than  the  echelon  flank. 

The  body  wave  phases  over  the  basin  are  composed  of  direct  and  multiplv 
reflected  waves,  head  waves,  and  diffracted  waves  from  the  corners.  At  close 
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velocity  seismograms  for  the  constant  velocity  basin  models  with  normal 
traction. 


Figure  2-9.  Ray  diagrams  for  the  piecewise  constant  velocity 

function  illustrating  the  corner  diffracted  arrivals 
excited  by  head  waves  on  the  interface  (critical 
angle  =  34°). 
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range  (out  to  3  or  4  km)  the  first  arrivals  are  direct  P-waves  and  are  the 
same  for  each  model.  Later  arrivals  are  masked  by  the  Rayleigh  wave.  At 
intermediate  range,  to  the  basin  edge,  the  arrivals  are  modified  by  the  flank. 
For  the  steep  flank  the  constructively  interfering  phases  are  virtually  iden¬ 
tical  to  the  layer  case  in  Fig.  2-7a.  However,  for  the  echelon  and  dipping 
flank, the  intermediate  range  arrivals  constructively  interfere  over  a  shorter 
segment  on  the  surface  but  yield  about  the  same  amplitude. 

Outside  the  basin  the  leading  body  wave  arrivals  are  due  entirely  to 
diffractions  from  the  corners  on  the  interface.  They  are  caused  by  the  head 
wave  traveling  along  the  interface  and  exciting  diffracted  waves  as  each 
corner  is  encountered.  Elementary  ray  diagrams  illustrating  the  process  are 

drawn  in  Fig.  2-9.  Further  magnification  of  the  seismograms  is  necessary 
to  clearly  identify  these  phases, as  they  are  weak  in  the  bedrock. 

The  later  phases  outside  the  basin  are  due  to  direct  body  and  surface 
waves  within  the  basin  interacting  with  the  corners.  The  strongest  is,  of 
course, the  transmitted  surface  wave,  but  it  is  preceeded  by  a  sequence  of 
weaker  Rayleigh  waves  excited  by  the  mode  conversion  of  direct  body,  head, 
and  diffracted  waves  at  the  surface  corner.  These  can  be  seen  as  leading 
phases  parallel  to  the  transmitted  Rayleigh  phase,  and  are  easily  identified 
by  viewing  obliquely  from  the  origin  along  this  strong  phase.  The  most  dis¬ 
tinct  example  is  for  the  steep  flank  in  Fig.  2-8,  where  direct  body  waves 
are  partially  converted  to  surface  waves  at  the  corner  starting  at  3  to  3.5 
sec.  Note  that  this  mode  conversion  also  excites  body  waves  outside  which 
are  too  weak  to  see  in  the  seismograms. 

2. 8  Basins  With  Linear  Velocity 

Synthetics  for  the  three  basin  models  with  piecewise  linear  velocity 
functions  and  normal  surface  traction  at  the  origin  are  pictured  in  Fig.  2-10. 
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In  these  cases,  the  seismograms  indicate  that  steep  and  echelon  faulted  flanks 
(Fig.  2-10a,b)  are  the  more  effective  reflectors  of  Rayleigh  waves,  while 
the  dipping  flank  (Fig.  2-10c)  is  a  poor  reflector.  In  terms  of  Rayleigh  wave 
transmission,  the  three  basin  flanks  appear  to  be  comparable.  These  conclusions 
are  consistent  with  the  previous  case  of  constant  velocity  basins. 

The  body  wave  phases  within  the  basins  are  similar  to  those  described  for 
the  layer  with  linear  velocity.  The  first  arrivals  out  to  about  4.5  km  are 
reflected  P-waves  in  the  basin.  From  4.5  km  to  the  basin  edge,  the  arrivals 
also  include  reflected  and  refracted  bedrock  waves,  as  well  as  diffracted  waves 
from  the  interface  edges.  However,  the  arrival  sequence  is  difficult  to  inter¬ 
pret  from  these  relatively  low  frequency  seismograms.  Better  resolution  would 
require  higher  frequencies,  hence  finer  grid  spacing  in  and  around  the  basin. 

In  both  the  echelon  faulted  and  dipping  flank  cases,  the  body  wave  amplitude 
decreases  beyond  5  km  due  to  scattering  and  mode  conversion  of  the  refracted 
P-wave  by  the  shallowing  interface.  The  predominant  body  wave  arrivals  follow¬ 
ing  the  P-waves  within  the  basin  are  refracted  S-waves.  These  are  clear  in  the 
steep  flank  case  (Fig.  2-10a),  but  are  complicated  by  mode  conversions  in  the 
other  basins.  In  the  magnified  seismograms  for  the  steep  and  echelon  faulted 
basins  (Fig.  2-10a,b),  the  refracted  S-waves  are  seen  to  reflect  very  weak 
Rayleigh  waves  back  into  the  basin  from  the  edge. 

The  body  wave  phases  outside  the  basins  are  due  to  bedrock  refractions, 
interface  edge  diffractions  and  mode  conversions  at  the  basin  edge.  The  first 
arrival  beyond  the  basins  is  the  refracted  P-wave.  The  refracted  S-wave  is 
the  next  significant  arrival  and  achieves  its  highest  amplitude  for  the  steep 
flanked  basin  and  is  less  pronounced  for  the  echelon  faulted  basin,  caused  in 
part  by  the  interference  of  edge  diffracted  waves.  For  the  dipping  flank  case, 
the  refracted  S-wave  is  very  weak,  probably  because  of  weaker  P-  to  S-mode 
conversion  on  the  sloping  interface. 
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The  surface  wave  phases  transmitted  outside  the  basin  are  caused  princi¬ 


pally  by  mode  conversions  at  the  basin  edges.  These  may  be  body  to  surface 
or  surface  to  surface  conversions.  The  strongest  surface  wave  is  the  direct 
conversion  (surface  to  surface),  as  was  the  case  for  the  piecewise  constant 
velocity  model.  In  the  present  case,  however,  there  are  significant  conver¬ 
sions  from  S-waves  to  surface  waves  at  the  basin  edge.  These  can  be  seen  in 
the  magnified  seismograms  as  the  weak  phases  preceding  the  direct  Rayleigh 
wave  phase  and  parallel  to  it.  They  are  most  coherent  for  the  dipping  flank 
(Fig.  2-10c),  probably  caused  by  multiple  reflection  and  conversions  of  the  body 
waves  as  they  reverberate  up  the  smooth  slope.  This  mechanism  is  described  in 
a  previous  study,  Wojcik  et  al.  (1981).  Conversions  from  P-waves  to  Rayleigh 
waves  are  not  indicated  by  the  seismograms.  Any  conversion  of  surface  waves  is 
obscured  by  the  transmitted  surface  wave  due  to  the  low  frequency  resolution 
of  the  model. 

2. 9  Basins  With  Tangential  Surface  Traction 

The  previous  seismograms  were  calculated  for  a  normal  surface  traction 
at  the  origin.  Synthetics  resulting  from  a  tangential  surface  traction  are 
illustrated  in  Fig.  2-11  for  the  three  basins  with  velocity  gradients 
described  above.  They  again  show  vertical  velocity  seismograms  normalized  by 
the  Rayleigh  wave  amplitude  300  m  from  the  basin  centerline  (second  seismogram). 
On  the  centerline  the  vertical  velocity  is  zero.  The  principal  Rayleigh  wave 
reflects  and  transmits  from  the  three  basin  edges,  as  described  above  in  Sec¬ 
tion  2.8,  but  the  body  waves  and  mode  converted  Rayleigh  waves  are  significantly 
different,  both  inside  and  outside. 

Within  the  basin,  P-waves  appear  as  before  but,  as  would  be  expected,  the 
S-wave  phase  is  clearly  stronger.  This  is  difficult  to  quantify  from  the 
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with  tangential  surface  traction. 


seismograms  but  is  evident  from  relative  amplitudes  of  the  converted  Rayleigh 
wave  (S-wave  to  Rayleigh  wave)  reflected  back  into  the  basin.  This  is  seen 
clearly  in  magnified  seismograms  for  the  steep  flanked  basin,  Fig.  2-lla,  and 
is  also  in  evidence  for  the  echelon  faulted  flank.  The  converted  Rayleigh 
wave  is  completely  absent  for  the  dipping  flank.  There  does  not  appear  to  be 
a  reflected  S-wave  in  the  basin  in  any  case,  probably  because  multiples  within 
the  basin  tend  to  obscure  weak  phases. 

The  principal  difference  between  normal  and  tangential  sources  appears  to 
be  body  wave  and  converted  phases  outside  the  basin.  For  the  steep  flanked 
basin.  Fig.  2-lla,  the  S-wave  from  the  interior  converts  to  a  Rayleigh  wave  and 
a  weaker,  somewhat  incoherent  S-wave  in  the  exterior.  In  contrast,  for  the 
echelon  faulted  and  dipping  flanks.  Fig.  2-llb,c,  the  interior  S-wave  is  trans¬ 
mitted  as  an  S-wave,  with  little  conversion  to  Rayleigh  waves.  These  trans¬ 
mitted  S-waves  are  virtually  nonexistent  in  the  case  of  normal  tractions.  Of 
course,  a  direct  comparison  of  amplitudes  between  normal  and  tangential  traction 
synthetics  is  not  valid  because  of  different  normalizations.  However,  relative 
scaling  of  S-wave  to  P-wave  amplitudes  in  each  case  confirm  that  the  tangential 
source  excites  much  higher  amplitude  S-waves  outside  the  basin  for  echelon 
faulted  and  dipping  flanks.  Reasons  for  these  differences  probably  involve 
diffractions  and  constructive  interference  of  multiple  arrivals  outside  the 
basin.  Detailed  analysis  is  beyond  the  scope  of  the  present  study. 

2. 10  The  Basin-Range-Basin  Model 

The  previous  finite  element  calculations  and  interpretations  serve  to 
characterize  a  variety  of  basin  interfaces,  velocity  functions  and  source 
types  for  a  single  symmetric  basin  in  a  half space.  In  this  section,  these 
results  are  extended  to  a  two-basin  model  with  an  intervening  1  km  mountain. 
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Vertical  velocity  synthetic  seismograms  are  drawn  in  Fig.  2-12  for  the  normal 


surface  traction  source.  The  wavelet  frequency  is  centered  at  1  Hz  in  these 
synthetics  to  minimize  data  storage  requirements  (rather  than  the  2  Hz  wavelet 
used  previously).  Amplitudes  are  normalized  by  the  Rayleigh  wave  maximum  on 
the  second  trace,  600  m  from  the  model  centerline  where  the  load  is  applied. 

Qualitative  features  of  the  seismograms  are  as  expected  from  the  previous 
basin  flank  synthetics.  The  principal  phase  is  the  Rayleigh  wave  which  is 
partially  reflected  and  transmitted  across  the  first  basin  edge.  The 
transmitted  surface  wave  travels  over  the  mountain  with  a  decreased  amplitude, 
due  to  scattering  losses  and  a  higher  characteristic  impedance,  and  is  trans¬ 
mitted  with  little  evidence  of  reflection  into  and  across  the  second  basin. 
Although  these  plane  strain  results  indicate  virtually  no  change  in  amplitude 
with  distance  (except  at  edges),  in  actuality,  for  a  point  source  the  Rayleigh 
wave  signal  would  decay  more  like  / 6/R  assuming  axisymmetric  geometry,  where  R 
is  distance  from  the  centerline  in  kilometers.  This  is,  of  course,  due  to  cylin¬ 
drical  divergence  of  the  surface  wave  (assuming  unit  amplitude  at  .6  km).  Thus, 
at  the  mountain  peak  (R  =  16  km),  the  amplitude  is  about  20  percent  of  that 
shown,  and  in  the  center  of  the  second  basin  (R  =  32  km),  it  is  reduced  to 
nearly  14  percent.  These  are  estimates  of  qualitative  behavior  only.  A  three- 
dimensional  calculation  is  necessary  to  quantify  the  true  amplitude  decay. 

Within  the  first  basin,  the  Rayleigh  wave  is  seen  to  undergo  multiple 
reverberations,  equivalent  to  reflections  from  the  opposite  basin  edge.  The 
second  reflection  at  23  seconds  and  8  km  in  the  magnified  synthetics  is  weak 
and  loses  coherence  in  the  background  noise;  however,  the  transmitted  wave  is 
still  significant.  Over  the  mountain,  the  direct  Rayleigh  wave  appears  to  be 
only  slightly  modified  by  the  topography  and  is  followed  15  seconds  later  by  the 
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much  weaker  transmitted  wave  from  the  reverberation  in  the  basin.  In  the  second 


basin,  the  direct  wave  transmits  across  the  left  edge  with  an  increased  ampli¬ 
tude  because  of  the  decrease  in  impedance.  There  is  no  coherent  indication 
of  direct  Rayleigh  wave  reflection  at  the  basin  edge.  This  is  probably  charac¬ 
teristic  of  incidence  from  a  high  impedance  to  low  impedance  medium,  but  also 
depends  on  details  of  the  edge  geometry.  Across  the  second  basin,  the  direct 
wave  propagates  uniformly. 

Body  wave  phases  over  the  left  basin  and  mountain  are  essentially  the  same 
as  those  found  in  the  steep  flanked  basin  calculation.  Fig.  2-10a,  with  allowance 
made  for  the  reduced  frequency  resolution  of  the  present  synthetics.  Arrivals 
on  the  mountain  consist  of  refracted  P-  and  S-waves,  a  mode  converted  Rayleigh 
wave  from  S-waves  in  the  left  basin,  and  the  direct  Rayleigh  wave.  In  contrast, 
over  the  second  basin  the  body  wave  arrivals  become  considerably  more  involved. 
Because  the  velocity  gradient  goes  to  zero  at  a  depth  of  5  km,  the  second  basin 
is  probably  in  a  shadow  for  refracted  P-  and  S-waves  through  the  bedrock.  There¬ 
fore,  the  first  body  wave  arrivals  are  shallower  refractions,  reflected  one  or 
more  times  from  the  free  surface.  From  the  magnified  synthetics,  the  first 
(refracted)  P-wave  is  seen  to  decay  very  rapidly  across  the  basin.  Subsequent 
body  wave  arrivals  are  much  stronger  and  appear  to  correlate  with  S-wave  and 
Rayleigh  wave  arrivals  at  the  left  edge.  That  is,  they  are  caused  by  diffrac¬ 
tion  and  mode  conversion  from  interaction  with  the  faulted  flank  on  the  left 
side  of  the  basin.  It  is  obvious  from  the  seismograms  that  to  understand  the 
complex  of  arrivals  in  the  second  basin  requires  both  higher  frequency  resolu¬ 
tion  and  a  detailed  ray  theory  analysis  to  discriminate  the  body  and  surface 
waves. 
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2.11  Selected  Particle  Motion  Plots 


To  facilitate  comparison  of  surface  wave  amplitudes  across  each  model, 
as  well  as  between  models,  particle  motion  plots  are  drawn  above  selected  seis¬ 
mograms  in  Figs.  2-13,14,15.  They  generally  show  retrograde  elliptical  orbits 
corresponding  to  the  Rayleigh  wave,  and  complicated  low  amplitude  motion  from 
the  body  waves. 

In  Fig.  2-13,  seismograms  and  orbits  are  shown  for  the  three  basin  types 
with  velocity  gradients  and  normal  traction  source  (described  in  Section  2.9). 
The  most  striking  feature  is  the  similarity  of  orbital  motion  between  each 
basin  type.  This  suggests  that  the  Rayleigh  wave  is  relatively  insensitive 
to  details  of  interface  geometry.  Wavelength  of  the  2  Hz  surface  waves  in  the 
basin  is  on  the  order  of  a  kilometer,  indicating  that  the  depth  of  significant 
motion  is  about  2  km.  Therefore,  in  the  deeper  basin  interior  away  from  the 
edges,  similar  orbits  would  be  expected.  However,  the  plots  show  that  this  is 
true  near  the  edges  as  well,  indicating  the  insensitivity  of  free  surface  Ray¬ 
leigh  waves  (in  contrast  to  ducted  Rayleigh  waves)  to  the  interface  shapes  and 
impedance  contrasts  assumed  here. 

Orbits  outside  the  three  basins  also  show  that  the  transmitted  surface 
wave  amplitude  is  relatively  independent  of  the  edge  geometry,  although  the 
details  of  motion  appear  to  differ  more  than  in  the  interior.  The  maximum 
amplitude  ratio  of  interior  to  exterior  waves  is  approximately  a  factor  of  three 
in  all  cases.  The  ellipticity  (ratio  of  vertical  to  horizontal  motion)  is  gen¬ 
erally  about  1.5.  Major  axes  of  the  elliptical  orbits  are  vertical,  except 
near  the  traction  source  and  very  close  to  the  basin  edge.  Near  the  source, 
the  orbit  is  skewed  because  strong  body  waves  are  superposed  on  the  surface 
wave.  However,  close  to  the  edge  the  skewness  is  probably  caused  by  the  local 
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c.  Dipping  flank 


Figure  2-13.  Particle  motion  plots  for  the  three  basin 
types  with  velocity  gradients  and  normal 
surface  traction. 
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change  in  wavespeed  in  the  horizontal  direction,  similar  to  the  effect  of 
anisotropy. 

Particle  motions  from  the  tangential  traction  source  are  plotted  for  the 
echelon  faulted  basin  in  Fig.  2-14.  They  show  a  more  erratic  motion  in  the 
first  two  orbital  plots  (compared  to  the  same  basin  in  Fig.  2-13)  because  of 
stronger  refracted  body  waves  from  the  tangential  source.  The  ellipticity  is 
found  to  be  in  the  range  1.4-1. 5  and  the  maximum  amplitude  ratio  of  interior  to 
exterior  waves  is  again  about  three  because  the  Rayleigh  wave,  once  generated, 
is  independent  of  source  type. 

In  Fig.  2-15,  particle  motion  is  plotted  for  the  basin -range-basin 
model.  Comparing  amplitude  in  the  left  basin  to  that  over  the  mountain  gives 
a  ratio  of  approximately  three,  although  there  appears  to  be  some  variability 
over  the  mountain.  The  amplitude  ratio  between  left  and  right  basins  is  about 
1.7.  As  for  any  of  these  plane  strain  models,  the  amplitude  ratios  should  be 
scaled  to  account  for  the  Rayleigh  wave's  cylindrical  divergence  from  a  more 
realistic  point  source. 
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SECTION  3 


THE  FULL  WAVE  FIELD  SOLVER 

This  section  describes  the  discrete  numerical  algorithm  used  to  solve 
the  time  domain  wave  field  equations  in  two-dimensional  inhomogeneous,  elastic 
and  acoustic  media.  The  physical  domain  is  a  rectangular  2-D  section  over 
which  the  Lagrangian,  small  displacement  computational  mesh  is  globally  carte¬ 
sian.  A  discrete  system  of  equations  is  derived  by  factoring  the  wave  field 
into  elastic  and  inertial  components  and  discretizing  each  separately.  A 
finite  element  displacement  approximation  over  rectangular  mesh  boxes  (ele¬ 
ments)  yields  elastic  force  resultants  at  the  mesh  points  (nodes),  while  a 
lumped  mass  approximation  at  the  mesh  points  provides  the  inertial  force 
resultants.  The  nodal  distribution  of  lumped  masses  (particles)  is  moved 
incrementally  in  time  according  to  the  discrete  impulse-momentum  equality, 
with  the  incremental  elastic  restoring  forces  at  each  node  accumulated  from 
contiguous  elements.  The  global  time  increment  is  chosen  smaller  than  the 
shortest  transit  time  across  any  element,  yielding  a  stable  explicit  integra¬ 
tion  algorithm  for  the  uncoupled  system  of  discrete  equations.  Large-scale 
modeling  applications  require  on  the  order  of  10^  nodes  entailing  approxi¬ 
mately  10*0  floating  point  operations  per  model.  Therefore,  coding  of  the 
algorithm  is  necessarily  designed  for  the  pipeline  and  parallel  capabilities 
of  existing  supercomputers,  specifically,  vectorization  on  the  CRAY-1. 

3 . 1  Background 

Solution  techniques  for  the  full  elastic/acoustic  wave  field  in  inhomo¬ 
geneous  media  generally  depend  on  some  form  of  discrete  algorithm.  The  usual 
basis  is  a  finite  difference  (FD)  or  finite  element  (FE)  spatial  discretiza¬ 
tion  over  an  appropriate  computational  mesh.  In  the  case  of  finite  differences, 
the  governing  partial  differential  equations  (PDE's)  are  typically  reduced 
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to  a  discrete  system  of  ordinary  differential  equations  (ODE's)  in  time  by 
means  of  difference  operators  for  spatial  derivatives  at  mesh  points.  In 
contrast,  with  finite  elements  the  particle  displacement  field  is  discretized 
via  element  shape  functions,  and  the  reduced  system  of  discrete  nodal  ODE's 
derived  by  enforcing  element  equilibrium.  Provided  the  differencing  operator 
and  element  shape  function  are  of  the  same  order,  both  FD  and  FE  discretiza¬ 
tions  result  in  similar  systems  to  be  solved,  e.g.,  Frazier  et  al.  (1973). 

Note  that  the  FE  method  does  not  explicitly  use  the  governing  PDE's  (as  does 
the  FD  method)  but  instead  satisfies  equilibrium  directly.  In  this  sense  it 
provides  a  more  fundamental  algorithm,  based  directly  on  elemental  mechanics 
rather  than  infinitesimal  mathematics  (which  are  abstracted  from  the  mech¬ 
anics).  This  is  best  appreciated  in  terms  of  ease  of  implementation  for 
complicated  media  and  shapes,  as  well  as  the  simplicity  of  elemental  solutions 
and  higher  order  generalizations. 

There  remains  time  integration  of  the  discrete  field  equations  by  means 
of  implicit  and/or  explicit  schemes,  whereby  the  ODE's  are  reduced  to  an 
algebraic  system.  These  equations  may  be  classified  as  hyperbolic  (wave-like) 
and/or  elliptic  (potential-like)  depending  on  the  time  scales  and  character¬ 
istic  lengths  in  the  problem.  The  terms,  hyperbolic-elliptic,  follow  from  the 
usual  classification  of  partial  differential  equations.  Briefly,  for  elastic/ 
acoustic  fields,  when  inertial  forces  dominate  elastic  forces,  the  system  is 
hyperbolic — but  when  elastic  forces  dominate,  it  is  elliptic.  In  a  hyperbolic 
system,  disturbances  at  one  point  do  not  influence  other  points  until  the  wave 
arrives.  However,  in  an  elliptic  system,  a  disturbance  at  any  point  influ¬ 
ences  every  other  point  without  time  delay.  Explicit  and  implicit  time  inte¬ 
grators  exhibit  similar  behavior.  Explicit  schemes  (hyperbolic  integrators) 
involve  the  independent  integration  of  each  nodal  equation  of  motion  by 
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virtue  of  the  hyperbolicity  of  the  problem.  That  is,  the  time  step  is  chosen 
small  enough  that  integration  (motion)  of  one  node  does  not  influence  an  adjacent 
node  during  one  time  step.  Implicit  schemes  (elliptic  integrators)  require 


the  simultaneous  integration  of  the  entire  system  of  equations,  assuming 
coupling  between  all  nodes  over  one  time  step,  i.e.,  elliptic  behavior.  This 
allows  the  use  of  a  somewhat  larger  time  step  than  for  the  explicit  case. 
Actually  all  nodes  are  not  coupled  in  hyperbolic  problems,  but  depend  on  each 
node's  sphere  of  influence,  proportional  to  the  time  step  chosen;  hence,  only 
those  within  the  sphere  are  coupled  in  an  implicit  scheme.  This  feature  makes 
an  implicit  solver  for  large  hyperbolic  systems  prohibitively  redundant. 

In  the  present  stud;.  ,  tine  domain  wave  propagation  of  pulses  is  the 
phenomenon  of  interest,  governed  by  hyperbolic  equations.  Because  the  require¬ 
ment  is  to  solve  very  large  systems  at  fast  cycling  rates  (mesh  or  element 
time  steps/second),  explicit  time  integration  of  the  discrete  equations  is 
essential.  The  choice  of  spatial  discretization,  finite  element  or  finite 
difference,  is  not  so  obvious  in  view  of  the  inherent  similarities  of  the 
resulting  systems.  However,  experience  with  the  implementation  of  both 
methods  shows  that  the  intuitive  and  mechanical  aspects  of  finite  element 
modeling  offer  clear  advantages.  Therefore,  an  explicit  finite  element 
algorithm  is  used  here,  and  is  described  in  the  remainder  of  this  section. 

3.2  Computational  Mesh 

The  first  consideration  in  implementing  a  large-scale  vectorized  wave 
solver  is  the  computational  mesh — i.e.,  the  distribution  of  discrete  field 
points  through  the  domain.  Given  a  rule  for  field  interpolation  between 
points,  the  mesh  must  be  capable  of  resolving  the  wave  pulse  in  space  and 
time,  as  well  as  inhomogeneities  in  the  medium.  In  this  study  it  must  also 
satisfy  a  large-scale  requirement — that  the  minimum  wavelength  resolved  be  a 
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small  fraction  (rather  than  a  multiple)  of  the  characteristic  dimension  of 
the  model.  Or  tne  computational  side,  the  mesh  must  possess  sufficient  regu¬ 
larity  to  permit  vectorization  of  the  algorithm  and  minimize  the  arithmetic 
and  storage  used  in  coordinate  processing.  Finally,  with  regard  to  physical 
modeling,  it  should  easily  represent  (from  a  user's  viewpoint)  complicated 
inhomogeneities  and  interfaces  to  a  degree  of  approximation  consistent  with 
its  frequency  resolution. 

The  simplest  computational  mesh  satisfying  these  requirements  is  carte¬ 
sian,  illustrated  in  Fig.  3-1  for  equidistant  grid  lines.  The  bold  contours 
are  a  fit  of  the  geologic  model  in  Fig.  2-1  showing  the  stepwise  resolution 
of  material  interfaces.  In  applications,  the  grid  density  necessary  for 
pulse  resolution  would  be  at  least  quadrupled,  with  sixteen  times  more  ele¬ 
ments — approximately  36,000  in  Fig.  3-1 — yielding  a  dense  matrix  of  (picture) 
elements  for  inhomogeneous  modeling.  This  identification  of  finite  elements 
with  picture  elements,  so-called  pixels  in  incremental  computer  graphics 
applications,  e.g. ,  Newman  and  Sproull  (1979),  provides  a  simple,  flexible 
scheme  for  constructing  large-scale  models  while  minimizing  user  and  computa¬ 
tional  overhead. 

To  propagate  a  wave  pulse  through  the  inhomogeneous  grid,  some  minimum 
number  of  elements  must  support  it,  depending  on  the  degree  of  field  inter¬ 
polation  function  used.  An  adequate  minimum  in  what  follows  is  ten  per  wave¬ 
length  of  the  highest  frequency  component  to  be  preserved  in  the  pulse. 
Therefore,  if  an  interface  in  the  cartesian  grid  is  composed  of  one  or  two 
element  steps  as  in  Fig.  3-1  (typically),  its  local  scale  is  equal  to  or  less 
than  one-fifth  of  the  shortest  wavelength.  A  somewhat  denser  support  of  the 
wave  over  the  mesh,  in  conjunction  with  inherent  smoothing  by  the  lumped  mass 
and  finite  element  discretizations,  makes  the  interface  practically  negligible 
in  terms  of  its  non-specular  scattering  effect  on  pulses. 
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Alternatively,  the  stepwise  interface  roughness  could  be  reduced  by 
using  conforming  elements,  e.g.,  quadralaterals  for  a  piecewise  linear  fit. 
However,  this  increases  arithmetic,  storage  and  the  level  of  element  approxi¬ 
mation,  as  well  as  complicating  the  vectorization  of  code  loops  necessary  for 
rapid  mesh  cycling.  For  these  reasons  the  simple  cartesian  mesh  with  stepwise 
resolution  of  inhomogeneities  will  be  used  exclusively  in  what  follows. 

3.3  The  Finite  Element  Approximation 

To  introduce  finite  element  discretizations  of  the  wave  field,  consider 
a  rectangular  domain  covered  by  a  cartesian  grid  of  rectangular  elements. 
Global  coordinates  are  (x,y)  and  local  (element)  coordinates  are  (£,r|)»  meas¬ 
ured  from  the  origin  at  (x^.y^)  of  the  ij  element.  These  are  illustrated  in 

T 

Fig.  3-2.  Assuming  that  the  material  displacement  vector,  d(x,y,t)  =  (u,v) 
is  regular  over  the  grid,  it  may  be  expanded  in  a  two-dimensional  time  depend¬ 
ent  Taylor  series  about  any  of  the  local  origins  as 


dO^+S,  yj+n,t)  =  d(x1,yj,t)  +  £dx  +  ndy^ 


i  i  i'j 


(3-1) 


where  0<£<Ax,  0<r|<Ay  and  derivatives  at  the  grid  points  are  denoted  by  sub¬ 
scripts,  e.g.,  d  .  When  the  grid  spacing  is  made  sufficiently  small  (but 
'xi 

finite),  only  certain  leading  terms  in  this  expansion  will  be  significant 
within  the  element.  If  these  principal  degrees  of  freedom  can  be  solved  in 
terms  of  nodal  displacements,  i.e.,  the  generalized  displacement  of  the  prob¬ 
lem,  then  local  continuum  analysis  may  proceed  on  this  basis.  Provided  the 
series  truncation  is  uniform  over  the  grid  of  local  expansions,  the  global 
field  is  continuous,  C^,  i.e.,  inter-element  displacements  are  compatible. 


-37- 


I 


The  simplest  local  interpolant  is  a  two-dimensional  version  of  the 
familiar  trapezoidal  rule,  which  follows  by  neglecting  curvature  terms,  i.e., 
second  and  higher  xi  and  y  derivatives  in  (3-1),  giving 


d(x  +-,y  +n,t)  -  d(x.,y.,t)  +  t,d  +  nd  +  5nd  (3-2) 

J  1  1  Xi  yj  xiyj 

For  this  truncation  of  the  Taylor  series  the  remaining  nodal  derivatives  are 
expressible  as  first  order  differences  of  the  element  nodal  displacements, 

=  (vvT  as 
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(3-3) 


where  n  runs  from  1  to  4  (clockwise  in  Fig.  3-2).  Introducing  the  nodal  dis¬ 
placement  vector, 

T  T 

6g  2  (u,v)  =  (u1,u2,u3,u4,v1,v2,v3,v4) 

the  element  expansion,  (3-2)  becomes. 


where 


and 
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(3-4) 


C  =  Ax-5  ,  n  =  Ay-n. 


Retention  of  higher  order  derivatives  in  the  truncation  requires  additional 
data  at  intermediate  points  on  or  in  the  element  for  a  unique  representation. 
In  view  of  the  minimal  arithmetic  involved  in  (3-4)  and  proven  utility  of  the 
trapezoidal  rule  in  numerical  analysis,  higher  order  elements  will  not  be 
considered  further. 
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With  an  approximation  (interpolant)  available  for  intra-element  displace¬ 
ments,  2-D  elasticity  theory  yields  elemental  stresses  and  strains.  The 
linear,  isotropic  stress-strain  relation  is 


cr=  (o 


a  )' 

nn’  Sn 


C  e. 


where  the  infinitesimal  strains  are 


(3-5) 


(£ 


KV  nn’ 


eCn> 


<ue*  V 


w 


C  is  the  constitutive  matrix  given  by 


(3-6) 


and  the  constants  are 


(3-7) 


Plane  strain:  ^  =  (3K+4G)/3  ,  c2  =  (3K-2G)/ (3K+4G)  ,  c3  =  G/^ 

Plane  stress:  ^  =  4G(3K+G)/ (3K+4G)  ,  c2  =  (3K-2G)/ (3K+G)/2  ,  c3  =  G/cj 

with  K  and  G  the  bulk  and  shear  moduli,  respectively .  Strains  are  evaluated 

T 

by  differentiating  d  *  (u,v)  in  (3-4)  giving 

e  *  B  6  (3-8 

-  e~e 

where  Bg  is  the  strain-displacement  matrix. 
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(3-9) 
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Therefore,  from  (3-5)  and  (3-8)  the  stress-displacement  relation  is 


a  =  CB  6 
e-e 


(3-10) 


The  final  step  is  to  calculate  the  generalized  elastic  and  inertial 
E  I 

forces,  f  =  f  +  f^,  i.e.,  the  equivalent  nodal  forces  corresponding  to  the 
generalized  displacements,  <5  .  This  is  accomplished  using  the  principle  of 
virtual  work.  For  a  virtual  displacement,  6^  of  the  nodes  the  internal  work 


is 


_  Ay  Ax  _ 

f  cj dA  =  /  /  6iB1CB  :  d£dn 


0  0 


~e  e  e-e 


and  the  external  work  of  the  nodal  elastic  forces  is  6  f  .  Equating  these 
expressions,  noting  that  5g  is  arbitrary,  yields 


f  =  K  5 
~e  e~e 


where 


Ay  Ax 

K,  =  /  /  BeCBedCdn 


(3-11) 


(3-12) 


0  0 

is  the  element  stiffness  matrix.  Similarly,  the  work  of  inertial  forces  due 
to  constant  mass  density,  0  is 

'T  ••  ~T  T 

l'd  (.'d)dA  =1  /  6  D  OD  6  d£dn 

A  '  '  0  0  'e  6  e~e 

^T  I 

and  that  of  the  corresponding  nodal  forces  is  -6  fg 


Equating  and  solving 


gives 


f  *  -M  6 
-e  e-e 


(3-13) 


where 


Ay  Ax 

M  =  ,-  /  /  De  Ded?dn 

0  0  e 


(3-14) 


is  the  element  mass  matrix. 
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The  symmetric  stiffness  and  mass  matrices  are  8x8  for  2-D  and  their 
closed  form  calculation  from  (3-12,  14)  is  straightforward,  involving  the 
integration  of  simple  polynomials  over  the  element.  The  resulting  stiffness 
matrix  is  fully  populated  while  the  mass  matrix  consists  of  two  4x4  blocks 
on  the  main  diagonal.  Both  are  given  in  Przemieniecki  (1968),  for  example, 
and  a  straightforward  numerical  evaluation  of  the  resulting  force  equations, 

(3-11,13)  requires  approximately  150  floating  point  operations  (adds,  multi- 

E  I 

plys,  etc.)  per  element  for  f  and  50  such  operations  for  fg. 

The  discrete  system  of  equations  governing  motion  of  the  single  finite 

£ 

element  follows  from  D'Alembert's  principle  (dynamic  equilibrium,  i.e.,  f  = 
!e  +  §e)  aS 


M  5  =  -K  6  +  g 

e-e  e-e 


(3-15) 


The  force  vector,  g^  represents  influence  of  adjacent  elements,  external 
forces,  etc.  In  order  to  solve  the  linear  system.it  is  very  useful  to  trans¬ 
form  the  equations  to  a  more  natural  displacement  basis  than  the  nodal  basis 
in  (3-15). 


3.4  The  Canonical  Finite  Element  Basis 

Element  stiffness  and  inertia  properties  can  be  defined  with  respect  to 
any  complete  set  of  displacement  modes.  These  may  be  nodal  displacements  as 
in  (3-11),  or  translations  and  some  combination  of  element  deformations,  e.g., 
extension,  shear,  dilatation,  rotation,  etc.  Although  the  cartesian  displace¬ 
ment  basis  used  above  is  ideal  for  derivation  of  the  governing  ODE's,  in 
terms  of  computational  efficiency  and  an  understanding  of  element  dynamics, 
it  is  not.  The  best  choice  is  the  canonical  finite  element  basis — defined 
here  to  be  a  displacement  basis  in  which  the  mass  matrix  is  diagonal  (by 
linear  transformation  rather  than  lumped  mass  approximation).  This  definition 
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requires  a  nontrivial  eigenvalue-eigenvector  analysis  of  (3-14)  to  determine 


the  basis.  Alternatively,  the  canonical  basis  may  be  defined  by  the  ortho¬ 
gonal  cartesian  displacement  modes  of  a  continuum  element:  translations 
Tx,Ty,  and  element  centered  extensions  Ex,Ey,  shears  Sx.Sy,  and  bends  Bx , By . 
These  are  illustrated  in  Fig.  3-3.  This  modal  definition  yields  a  convenient 
mechanical  interpretation. 

The  canonical  basis  is  found  most  directly  by  inspection  of  the  modes 
pictured  in  Fig.  3-3.  Writing  the  modal  vector  as 

6^  =  (Tx,Ex,Sx,Bx,Ty,Ey,Sy,By)T  (3-16) 


the  normalized  orthogonal  transformation,  S  between  nodal  and  modal  basis 
vectors  is 


5  =  So 

~e  ~m 

where  S  is  unitary 


(s~1=sT) 


1-1-1  1 
1  1-1-1 
1111 
1-1  1-1 


and  its  nonzero  blocks  are 


(3-17) 


(3-18) 


The  columns  of  this  transformation  constitute  the  modes  (eigenvectors),  i.e., 
the  normalized  nodal  displacements  ordered  clockwise  from  the  origin  in  Fig. 
3-3.  In  the  new  basis,  element  displacement  and  stress-strain  relations  are 


d  =  D  5  =  D  5  =DS6  ,  (3-19) 

-  e~e  m~m  e  ~m 

e=B6=B6=BS6  ,  (3-20) 

~  e~e  m~m  e  ~m 

a  -  Cc  =  DB  S6  .  (3-21) 

~  ~  e  ~m 

Evaluating  the  modal  strain-  and  stress-displacement  matrices,  the  strain  and 
stress  equations  are 
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(3-22) 
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(3-23) 


which  demonstrates  that  translation  (columns  1  and  5)  makes  no  contribution 
to  the  stress-strain  state,  while  extension  (columns  2  and  6)  and  shear 
(3  and  7)  produce  uniform  states  over  the  element,  and  bending  (4  and  8) 
adds  linear  variation. 

r  T 

To  express  the  stiffness  and  inertia  matrices,  K  =  J  B  CB  dA  and 

m  ,  m  m 

T  A 

M  =  p  /  D  D  dA  respectively;  in  the  new  basis  it  is  only  necessary  to  set 
in  .mm 

A 

B  =  B  S  and  D  =  D  S  and  evaluate  the  integrals,  giving 
m  e  m  e 


K  =  STK  S 
m  e 


0  0  0  0 

aL  0  0  0 

0  a^c3  0  0 

0  0  bx  0 

0  0  0  0 

c2  0  0  0 

0  c3  0  0 

0  0  0  0 


(3-24) 
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(3-25) 


where  m  =  pAxAy  is  the  element  mass  and  =  Ay/Ax,  a^  =  Ax/Ay,  =  (a^+ 
ajC ^)/3,  b2  =  (a2+a^c^)/3.  This  demonstrates  that  the  canonical  basis  does 
indeed  diagonalize  the  mass  matrix,  as  well  as  significantly  reducing  the 
number  of  nonzero  entries  in  the  stiffness  matrix  (from  69  to  10)  due  to  the 
simple  modal  stress  states. 


3.5  Canonical  Element  Dynamics 

In  the  canonical  basis  the  ODE's  governing  element  dynamics  become 


M  6  •  -K  J  +  g 

m~m  m~m  xm 


(3-26) 


With  diagonal  and  Km  sparse,  this  representation  allows  closed  form  solu¬ 
tions.  In  particular,  the  homogeneous  system  (setting  g^  =  0)  reduces  to  the 
following  set  of  equations  for  the  element  centered  free-vibration  modes. 


(3-27) 
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where  only  the  bending  modes  are  uncoupled.  Assuming  time  harmonic  solutions, 
i.e.,  ae^ut  and  solving  for  the  unknown  frequencies  and  modal  coupling, 
Extensional  vibrations 


uEx  =  “Ey  *  2  y<3c1/2m)(a1+a2^(a1-a2)  +4c2' 
Ey/Ex  *  (a2~ai  (a\~a2^  (2c2^ 


Shearing  vibrations 


(3-28) 


0  ;  Sy/Sx  =  -a. 


S“  '  I  x  v  (3G/m) (a,+a  )  ;  Sy/Sx  =  a 


(3-29) 


Bending  vibrations 


x  *  c . b , /m  , 


-  xc,b-/ra 


wBx  2  *  ^By  2 


(3-30) 


Specializing  these  solutions  to  a  square  element  (Ax=Av,  hence  a^=a2=l, 
b^*b2* (l+c^) /3)  yields  the  simple  picture  of  free  element  dynamics  shown  in 
Fig.  3-4.  The  zero  frequency  shearing  vibration,  Sy/Sx  ■  -1,  is  seen  to 
yield  rigid  body  rotation,  hence  there  are  only  five  deforraational  modes. 

The  pure  shear  mode,  Sy/Sx  *  1  and  pull-push  extensional  mode,  Ey/Ex  =  -1 
have  the  lowest  natural  frequency,  w^=*3/2  cg/Ax;  while  the  bending  modes, 
Bx  and  By  and  the  pull-pull  extensional  mode,  Ey/Ex  =  1  have  the  highest, 


-jgx  =*  *'(kz+l)/2  uj^  and  ^Ex  =  »kz-l  ^  ta£x  ><  Wgx  when  k  ><  v2),  where  k  = 


*  ^  2  —  1  uJ,,  (out 


0^/0^  is  the  wavespeed  ratio.  When  c^  =  0  (k-*°°)  in  an  acoustic  medium  these 


natural  frequencies  reduce  to  ^  =  0,  =  ✓  3/2  c  /Ax  and 

N  BX  P 


y'3/2  c  /  Ax . 
P 


The  natural  period,  T^  =  2^/^  may  be  written  as  T^  =  5.13k  Ax/c^.  In  geo¬ 
logic  media  for  example,  where  k  =  /3  is  a  good  approximation,  the  natural 
element  period  is  T^  =  8.89Atg,  with  Atg  denoting  p-wave  transit  time  across 
the  element.  Thus,  in  applications  is  a  clear  upper  bound  to  the  frequency 
resolution  of  a  finite  element  mesh  based  on  this  element.  In  operation,  if 


the  mesh  is  driven  externally  at  frequencies  approaching  or  exceeding 


r 


individual  elements  will  resonate  at  their  modal  frequencies,  introducing  a 
spurious  high  frequency  ringing  to  the  wave  field.  This  is.  of  course,  equiva¬ 
lent  to  the  numerical  dispersion  generally  associated  with  discrete  mathema¬ 
tical  systems  or  the  ringing  of  any  mechanical  (force-motion)  transducer. 

3.6  Explicit  Integration  and  the  Lumped  Mass  Approximation 

■Jith  the  element  ODE's  known,  the  equations  must  be  assembled  (coupled) 
and  solved  by  numerical  integration  in  time,  specifically  by  advancing  the 
global  system  incrementally,  at  timestep  At.  However,  the  degree  of  system 
coupling  depends  on  the  extent  of  dynamic  interaction  between  nodes  during 
a  single  timestep.  This  is  measured  by  the  radius  of  a  node's  circle  of 
influence  (sonic  circle),  r  =  c^t  shown  graphically  in  Fig.  3-5.  If  It  is 
less  than  the  shortest  element  transit  time,  i.e.,  if  the  circle  is  contained 
within  a  quad  of  elements  (the  smallest  circle  in  Fig.  3-5),  then  no  global 
interaction  occurs  because  all  nodal  motions  are  uncoupled  during  one  time- 
step.  Such  behavior  follows  from  hyperbolicitv  of  the  system  and  provides 
the  basis  for  the  explicit  integration  scheme  used  here.  For  larger  It  the 
equations  for  elements  within  and  including  the  sonic  circle  must  be  assembled 
and  solved  simultaneously  using  an  implicit  integration  approach. 

To  proceed, it  is  useful  to  write  the  assembled  system  of  model  equations 
as 

MU  =  -KU  +  G  (3-31) 

where  U  is  the  global  nodal  displacement  vector,  M  and  K  are  the  assembled 
mass  and  stiffness  matrices,  and  G  is  the  global  external  force  vector.  The 
symmetric  banded  structure  of  M  and  K  depends  on  the  nodal  ordering  of  U. 

They  are  assembled  by  first  switching  rows  and  columns  in  the  element 
matrices  and  inflating  with  rows  and  columns  of  zeros  so  the  ordering  and 
spacing  of  element  nodal  displacements  correspond  to  the  global  arrangement 
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Figure  3-5.  A  schematic  of  dynamic  interaction  between  nodes 
for  various  it  measured  by  the  central  node's 
circle  of  influence  (sonic  circle). 


in  U;  and  second,  summing  the  inflated  mass  and  stiffness  matrices  to  give 
M  and  K,  respectively. 

For  explicit  integration  the  internal  force  vector,  F  =  -KU  need  not  be 
assembled  but  is  instead  determined  by  accumulating  nodal  forces  via  an 
element  by  element  evaluation  of  the  element  stiffness  equation,  (3-11).  With 
the  right  hand  side  of  (3-31)  known,  the  incremental  form  of  the  equations 
is 

M  ^  =  F  +  G  ,  (3-32) 

whence  the  equations  for  velocity  increments  over  the  timestep  become 

MAU  =  At (F+G)  (3-33) 

This  is  recognized  as  the  incremental  impulse-momemtum  equality.  The  result¬ 
ing  integration  algorithm  includes  the  following  steps, 

1.  AU  =  M-1At (F+G) 

2.  U  =  U  +  AU 

3.  A U  =  At  U 

4.  AF  =  -KAU 

5.  F  =  F  +  AF 

6.  G  =  G  +  AG  (3-34) 

where  the  external  force  increments,  AG,  are  assumed  given  in  the  problem 
definition.  Returning  to  step  1,  the  process  is  repeated  for  the  next  and 
subsequent  timesteps.  In  this  algorithm  nodal  velocities  should  be  con¬ 
sidered  mean  values,  hence  are  formally  evaluated  at  the  timestep  midpoint, 
while  forces  and  displacements  are  evaluated  at  each  timestep.  However, 

in  applications  to  linear  problems  no  distinction  need  be  made.  Note  that 
nodal  displacements  are  not  calculated  during  the  integration,  as  only 
increments  are  used  to  find  AF  thus  F. 
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To  solve  for  AU  in  (3-33)  it  is  necessary  to  premultiply  by  the  mass 

matrix  inverse,  which  requires  assembly  and  inversion  of  the  global  mass 

matrix  and  complicates  the  algorithm  considerably.  The  explicit  scheme  averts 

these  global  matrix  operations  by  diagonalizing  the  mass  matrix  using  the 

lumped  mass  approximation.  On  the  element  level  this  amounts  to  replacing 

M  defined  in  (3-14)  by  ml/4  where  I  is  the  identity  matrix  and  m  is  the 
e 

element  mass,  pAxAy.  Assembling  the  global  mass  matrix  yields  a  diagonal 
form,  with  each  nodal  entry  having  1/4  the  mass  of  contiguous  elements. 

From  the  previous  analysis  of  element  dynamics  the  effect  of  mass  lump¬ 
ing  is  easily  found.  The  relation  between  modal  and  nodal  element  mass 
T  T 

matrices  is  M  =  S  M  S  so  that  (M  ),  .  =  S  (ml/4)S  =  ml/4,  whence, 

m  e  m  lumped 

(M  ),  ,  =  (M  ).  ,.  Comparing  M  in  (3-25)  to  (M  ),  shows  that 

m  lumped  a  lumped  m  m  lumped 

the  lumped  mass  is  3/4  of  the  consistent  mass  for  extensional  and  shear 
modes  and  9/4  for  bending  modes.  Substituting  into  (3-28)  -  (3-30)  gives 
(ui) lumped^  =  ~  f°r  extension  and  shear,  and  “  2/3  - 

0.667  for  bending. 

3. 7  Further  Development 

This  final  form  of  the  explicit  finite  element  algorithm  could  have  been 
written  more  directly  than  the  above  derivation  indicates,  with  only  minor 
recourse  to  theory  (or  a  textbook)  for  the  stiffness  matrix.  However,  the 
requirements  of  this  study — to  solve  very  large  problems  as  efficiently  as 
possible  on  supercomputers — warrant  a  closer  look  at  the  basic  algorithm.  The 
above  exposition  and  theory  are  the  basis  for  a  more  effective  formulation  of 
cartesian  finite  elements. 

The  theory  is  developed  here  to  a  level  commensurate  with  the  algorithm 
embodied  in  the  FORTRAN  code  applied  in  Section  2.  However,  the  theory  can  be 


carried  further.  In  particular,  the  idea  of  applying  a  linear  transformation 


to  diagonalize  the  element  mass  matrix  can  be  applied  to  the  element  stiffness 
matrix  yielding  the  free  vibration  modes  in  Fig.  3-4,  which  further  reduce  the 
operations  count.  This  provides  a  bonus  in  that  mathematically  consistent  fre¬ 
quency  independent  damping  can  now  be  selectively  defined  for  each  of  the  ele¬ 
ment  vibration  modes,  e.g.,  the  pure  dilatational  and  pure  shear  modes.  This 
is  equivalent  to  Rayleigh  damping  in  the  canonical  stiffness  basis  and  appears 
to  be  a  new  implementation  of  damping  in  explicit  finite  element  codes  (for 
which  frequency  independent  damping  has  always  been  a  problem).  Furthermore, 
the  global  assembly  of  element  equations  in  cartesian  coordinates  is  facilitated 
by  linear  transformation  to  diagonal  mass  and  stiffness  matrices  because  the 
inflated  transformation  matrices  have  very  simple  inverses. 

These  concepts  will  be  explored  fully  in  the  second  part  of  this  research 


study  involving  the  effects  of  material  damping  and  three-dimensional  modeling. 
Results  will  be  reported  in  a  second  and  final  report. 


SECTION  4 


SUMMARY  AND  CONCLUSIONS 

This  study  was  motivated  by  the  need  to  understand  free-field  ground 
motions  over  laterally  inhomogeneous  models  of  the  Earth's  crust — specifically, 
graben  basins  typical  of  the  Basin  and  Range  province.  To  address  this  ques¬ 
tion  on  the  physical  scale  of  a  nominal  graben  basin,  i.e.,  16  km  long,  2  km 
deep  with  wavelengths  as  small  as  500  meters,  requires  a  level  of  wave  propaga¬ 
tion  computing  power  approximately  one  order  of  magnitude  greater  than  scalar 
computers  allow  (e.g.,  CDC  7600).  This  scale  of  computation  necessitates  using 
a  powerful  vector  computer  on  the  order  of  a  CRAY-1  or  CYBER  205  to  be  economi¬ 
cally  feasible.  However,  in  addition  to  a  supercomputer,  an  efficient  vector¬ 
ized  algorithm  is  necessary  to  approach  the  theoretical  limit  of  performance. 
This  report  has  addressed  both  the  seismological  calculations  in  three  basin 
configurations,  executed  on  the  CRAY-1,  and  the  finite  element  theory  used  to 
formulate  them. 

4. 1  Summary  of  Calculations 

Synthetic  seismograms  over  the  basin  flank  models  are  characterized  by 
the  relative  strength  or  weakness  of  a  few  principal  phases.  Typical  travel 
time  curves  for  these  phases,  obtained  from  the  synthetic  seismograms,  are 
drawn  in  Fig.  4-1.  Inside  the  basin,  the  first  coherent  arrivals  are  P-waves, 
designated  Ph  1  in  the  figure,  including  direct  refracted  as  well  as  ref lected  and 
diffracted  waves  from  the  interface.  The  second  arrivals,  Ph  2,  are  S-waves 
following  similar  travel  paths  within  the  basin.  These  S-waves  may  have  suf¬ 
ficient  amplitude  to  reflect  significant  Rayleigh  waves,  Ph  3,  by  mode  conver¬ 
sion  at  the  basin  edge.  The  direct  Rayleigh  wave,  Ph  4,  is  the  strongest  phase 
in  all  cases  and  reflects  a  weaker  Rayleigh  wave,  Ph  5,  back  into  the  basin. 
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synthetic  seismograms. 


Outside  the  basin,  the  principal  body  wave  phases  include  a  refracted  P-wave, 

Ph  6,  and  a  refracted  S-wave,  Ph  7,  probably  including  a  transmitted  S-wave 
from  the  interior.  The  next  arrival  is  a  mode  converted  Rayleigh  wave,  Ph  8, 
from  the  S-wave  within  the  basin.  This  phase  is  typically  offset  in  time, 
probably  due  to  a  strong  phase  shift  in  the  conversion  process  near  the  edge. 
The  strongest  phase  outside  the  basin  is  the  transmitted  Rayleigh  wave,  Ph  9, 
from  the  direct  wave  in  the  basin. 

The  phases  shown  in  Fig.  4-’  are  incomplete  in  that  there  are  a  number  of 
other  very  weak  arrivals  discernable  in  the  seismograms,  but  difficult  to 
identify  visually  as  to  source  and  type.  Also,  the  actual  ray  paths  and  mode 
conversion  mechanisms,  i.e.,  multiple  reflections,  diffractions,  etc.,  of  the 
identified  phases  require  more  careful  study  to  quantify.  Therefore,  the  above 
descriptions  should  be  considered  very  qualitative,  subject  to  a  careful  ray 
analysis  to  identify  single  and  multiple  paths,  and  higher  frequency  finite 
element  analysis  to  allow  better  separation  of  phases  in  time  for  amplitude 
estimates. 

The  basin-range-basin  model  and  calculations  extend  the  physical  scale 
of  modeling  to  include  the  interaction  of  two  basins  separated  by  a  mountain. 
This  demonstrates  the  practicality  of  very  large-scale  calculations,  i.e., 
greater  than  100,000  elements,  on  the  CRAY-1.  The  resulting  synthetic  seismo¬ 
grams  show  phases  very  similar  to  those  described  above  for  the  basin  flanks 
and  do  not  indicate  any  strong  modification  by  the  mountain.  However,  in  the 
second  basin  the  phases  are  quite  complex,  caused  by  incident  body  and  surface 
waves  interacting  with  the  basin's  leading  edge.  This  is  in  contrast  to  the 
simple  surface  source  applied  to  the  left  basin  centerline  and  requires  addi¬ 
tional  analysis  to  quantify  the  resulting  phases.  The  strongest  phase  is 
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certainly  the  direct  Rayleigh  wave  which  does  not  appear  to  be  strongly  affected 
by  the  mountain  and  couples  efficiently  into  the  second  basin.  Mode  converted 
S-waves  (Rayleigh  to  S-wave)  from  the  leading  edge  are  the  strongest  body  wave 
phases  within  the  second  basin. 

The  effects  of  basin  interface  geometry,  source  type,  and  velocity  functions 
can  be  briefly  summarized.  The  principal  observation  is  that  basin  edges  are 
efficient  transmitters  of  Rayleigh  waves,  apparently  independent  of  the  edge 
details.  Rayleigh  wave  reflection  depends  strongly  on  the  edge  slope,  with 
steep  edges  better  reflectors  than  dipping  edges.  Body  waves  inside 
a  basin  are  affected  by  interface  shape  to  some  extent,  although  less  than  might 
be  expected.  Higher  frequency  analysis  is  needed  to  clarify  details.  Shear 
waves  outside  a  basin  depend  markedly  on  the  edge  geometry  due  to  mode  conver¬ 
sion  of  waves  in  the  interior.  These  conversions  favor  S-waves  for  the  steep 
flank  and  Rayleigh  waves  for  the  dipping  flank.  The  above  observations  apply  to 
the  normal  surface  source.  The  tangential  surface  source  generally  excites 
higher  amplitude  shear  waves  within  the  basin  favoring  stronger  mode  conver¬ 
sions  (S  to  S  and  S  to  Rayleigh)  at  the  basin  edge.  Comparing  the  two  load¬ 
ing  cases,  the  obvious  feature  outside  the  basins  is  the  strong  S-wave  for 
the  tangential  source  in  the  echelon  faulted  and  dipping  flank,  with  little 
or  no  converted  Rayleigh  wave.  This  situation  is  reversed  for  the  steep 
flanked  basin.  In  terms  of  velocity  functions,  synthetics  for  the  constant 
velocity  cases  are  qualitatively  similar  to  the  velocity  gradient  cases  but 
lose  much  of  the  detail.  They  are,  therefore,  instructive  but  not  quantita¬ 
tively  useful  in  modeling  real  structures. 

4. 2  Summary  of  Finite  Element  Theory 

Solving  for  seismic  waves  in  completely  inhomogeneous  elastic  media 
requires  some  form  of  numerical  analysis,  based  either  on  ray  or  wave  theory. 
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However,  only  wave  theory  provides  complete  solutions  (in  the  sense  that  all 
wave  types  are  included) .  The  most  common  implementations  of  wave  theory 
use  finite  difference  or  finite  element  discretizations.  For  reasons  of  effi¬ 
ciency  and  accuracy  (because  they  provide  a  closer  link  to  the  fundamental 
mechanics  of  the  seismic  medium),  finite  element  methods  are  preferred  here. 

Time  integration  of  the  finite  element  equations  is  performed  using  an  explicit 
approach,  which  is  very  efficient  for  transient  phenomena,  in  contrast  to 
implicit  methods.  The  two-dimensional  computational  mesh  is  a  simple  carte¬ 
sian  grid  with  a  stepwise  resolution  of  material  interfaces. 

The  finite  element  algorithm  is  implemented  in  FLEX,  a  standard  FORTRAN 
computer  code  designed  for  highly  efficient  processing  on  scalar,  as  well  as 
vector  computers.  A  brief  summary  of  the  algorithm  follows.  A  finite  element 
displacement  approximation  over  rectangular  elements  yields  elastic  force  resul¬ 
tants  at  the  nodes,  while  a  lumped  mass  approximation  provides  the  inertial 
force  resultants.  The  nodal  distribution  of  lumped  masses  is  moved  incrementally 
in  time  according  to  the  discrete  impulse-momentum  equality,  with  the  incremental 
elastic  restoring  forces  at  each  node  accumulated  from  contiguous  elements.  The 
global  time  increment  is  chosen  smaller  than  the  shortest  transit  time  across 
any  element,  yielding  a  stable  explicit  integration  algorithm  for  the  uncoupled 
system  of  discrete  equations.  A  normal  impedance  type  boundary  condition  is 
included  to  mimic  a  transmitting  boundary  (approximately). 

The  theory  described  in  Section  3  is  aimed  at  reducing  the  number  of  arith¬ 
metic  operations  required  to  evaluate  the  finite  element  equations.  The  standard 
engineering  implementation  of  finite  elements  allows  general  skewed  elements,  but 
when  rectangular  elements  only  are  used,  much  simpler  closed  form  expressions 
result.  In  particular,  it  is  found  that  linear  transformations  exist  which 
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diagonalize  the  element  mass  and  stiffness  matrices,  thereby  reducing  substanti¬ 
ally  the  operations  count.  These  transformations  are  equivalent  to  well-known 
deformational  modes  in  elasticity  theory.  An  added  benefit  of  the  diagonalizing 
transformation  on  the  stiffness  matrix  is  that  frequency  independent  damping  may 
be  applied  to  the  pure  dilatational  and  pure  shear  deformational  modes.  This  is 
accomplished  in  a  mathematically  consistent  fashion  and  adds  little  computa¬ 
tional  overhead. 

4. 3  Conclusions 

The  conclusions  to  be  drawn  from  this  one-year  research  effort  concern: 
large-scale  seismic  modeling  on  supercomputers;  graben  basin  response  in  the 
Basin  and  Range  province;  and  explicit  finite  element  algorithms.  The  princi¬ 
pal  conclusion  is  that  complete,  linear  elastic,  high  frequency  wave  fields  in 
inhomogeneous  2-D  media  can  be  economically  calculated  on  a  large  physical 
scale  using  the  CRAY-1  supercomputer  and  a  vectorized  elastic  wave  solver. 

For  example,  comparisons  made  elsewhere  between  this  approach  and  more  conven¬ 
tional  wave  solvers  on  the  CDC  7600  show  a  speed  increase  of  100  or  more  and  a 
size  increase  of  at  least  10.  Such  a  jump  in  computational  power  allows  the 
analyst  to  solve  problems  which  could  not  be  contemplated  before,  and  at  rela¬ 
tively  low  cost. 

The  focus  of  the  computational  effort  was  seismic  resoonse  of  three  basins 
for  normal  and  tangential  surface  sources.  The  suite  of  calculations  showed 
a  number  of  differences  between  a  steep  faulted  flank,  an  echelon  faulted  flank 
and  a  dipping  flank.  They  also  showed  that  full  field  synthetics  can  be  nearly 
as  complicated  as  real  data  to  interpret,  requiring  interpretational  tools  like 
ray  tracing  to  identify  phases.  Synthetics  were  generated  by  convolving  a  2  Hz 
wavelet  with  a  Green's  function  from  the  finite  element  run,  but  no  comparisons 
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with  other  time  histories  were  made.  Only  a  single,  multilinear  velocity  func¬ 
tion  was  used  in  the  modeling  (besides  the  piecewise  constant  function  which  was 
too  idealized)  and  it  would  be  useful  to  examine  variations  on  the  velocity 
function  for  a  single  basin  model.  No  attempt  was  made  to  include  material 
attenuation,  but  the  finite  element  implementation  of  damping  was  developed 
and  will  be  implemented  in  the  second  phase  of  this  effort.  The  normal  impedance 
boundary  absorber  with  element  grading  near  the  boundary  appeared  to  suppress 
boundary  reflection  for  the  cases  considered.  However,  a  better  theory  would 
help  to  reduce  problem  size,  especially  in  3-D  applications  to  be  examined  in 
the  second  phase. 

The  major  result  of  the  finite  element  theoretical  work  was  that,  in  the 
context  of  cartesian  finite  element  grids,  the  governing  equations  can  be  greatly 
simplified,  and  formulated  so  as  to  admit  selective  seismic  damping  on  pure  shear 
and  dilatational  modes  of  deformation.  In  addition,  the  application  of  linear 
transformations  on  the  global,  as  well  as  local,  level  offer  the  possibility 
of  further  reducing  operation  counts  and  increasing  execution  speed  of  the  fi¬ 
nite  element  algorithm.  This  will  be  examined  in  the  second  phase  for  3-D  appli¬ 
cations. 
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